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Abstract. Mixed positive and negative feedback loops are often found in biological 
systems which support oscillations. In this work we consider a prototype of such 
systems, which has been recently found at the core of many genetic circuits showing 
oscillatory behaviour. Our model consists of two interacting species A and B, where A 
activates not only its own production, but also that of its repressor B. While the self- 
activation of A leads already to a bistable unit, the coupling with a negative feedback 
loop via B makes the unit frustrated. In the deterministic limit of infinitely many 
molecules, such a bistable frustrated unit is known to show excitable and oscillatory 
dynamics, depending on the maximum production rate of A which acts as a control 
parameter. We study this model in its fully stochastic version and we find oscillations 
even for parameters which in the deterministic limit are deeply in the fixed-point 
regime. The deeper we go into this regime, the more irregular these oscillations 
are, becoming finally random excitations whenever fluctuations allow the system to 
overcome the barrier for a large excursion in phase space. The fluctuations can no 
longer be fully treated as a perturbation. The smaller the system size (the number 
of molecules), the more frequent are these excitations. Therefore, stochasticity caused 
by demographic noise makes this unit even more flexible with respect to its oscillatory 
behaviour. 



PACS numbers: 87.10.Mn,87.16.dj,87.16.Yc,87.18.Cf 
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1. Introduction 

Biological systems such as populations of interacting organisms or genetic networks 
are inherently non-linear due to feedback loops and stochastic due to a finite number 
of contributing ingredients and a discrete set of reaction events. One of the various 
origins for stochastic behaviour is demographic noise caused by fluctuations in the 
population size. Although often ignored in the past, it is now well established that 
such inherent stochasticity is not just a small correction to deterministic, infinite- 
population solutions, but it may lead to a number of qualitatively new effects. As 
already demonstrated, demographic noise can lead to temporal oscillations for predator- 
prey models [TJ, counter-intuitive reversions in the evolutionary process for replicator 
dynamics [2], or persistent spatial patterns [3]. In gene expression models, it has been 
shown for a single self-regulating gene that the solution of the stationary state reveals 
qualitative deviations from the deterministic solution like anticooperative behaviour in 
the limit of slow binding and unbinding rates jU [5] . In networks of mutually interacting 
genes, extinction and resurrection events for small population size can lead to additional 
fixed points in the attractor landscape [6]. In the repressilator model, a prototype of 
genetic oscillators, coherence resonance has been observed due to stochastic fluctuations 
[7] , and a modification of frequency, amplitude and parameter regime where oscillations 
occur has been demonstrated [8]. Similar effects of intrinsic noise have been seen in 
stochastic delay systems of genes [9]. 

The aim of this work is to analyse the effect of finite-population demographic noise 
on a so-called bistable frustrated unit (BFU). This model consists of two interacting 
species A,B. Species A activates its own production and without coupling to the second 
species B, the self-loop of A would correspond to a bistable system. But A here also 
activates the production of its repressor B. The unit of the two loops was then termed 
frustrated in analogy to antiferromagnetic frustrated couplings pU]. The reason is that 
A gets conflicting input, activation from itself and repression via the activation of the 
second species B. Throughout this paper we use "bistable frustrated unit" as a name for 
our model. What makes this simple model particularly interesting is that by tuning the 
ratio of half life of the two species, one can obtain fast dynamics for the activating species 
A and slow dynamics for the repressing species B, similarly to the fast and slow degrees 
of freedom of FitzHugh-Nagumo units as models for neural networks [HI H2]- This 
defines an intrinsic ratio of time scales with implications for the shape of amplitudes, 
the form of limit cycles, and the probability distribution in phase space. 

The deterministic (infinite-population) version of this model was originally proposed 
in Ref. [10] as an effective, coarse-grained description of genetic circuits in which a 
system, bistable due to a positive feedback loop, is coupled to a negative feedback loop 
which adds frustration to the system. The model has been shown to have a variety of 
oscillatory behaviour, which is easily tunable by the model parameters. In addition to 
the limit-cycle oscillatory regime, two fixed-point regimes have been recently identified 
in this model |13j . 
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In this paper we shall analyse a fully stochastic version of this model with a finite 
number of particles of each species. We shall see that, in contrast to the deterministic 
limit, cycles (, i.e., oscillations) appear even if the parameters of the model are deeply in 
the deterministic fixed-point regime. The existence of such "quasi- cycles" , as we shall 
call them in order to stress their difference to usual, deterministic cycles, indicates that 
the range of parameters for which the stochastic system exhibits oscillatory behaviour 
is much broader, as argued in Ref. [H] for a similar system. This suggests that fine- 
tuning of the parameters is less critical. This observation is relevant for real-world 
genetic circuits when they are realized as single units and exhibit oscillations such as 
the circadian clock [15], because it means that a finite number of particles can make 
them more robust to damage, e.g., caused by mutations changing protein production 
rates, binding affinities etc. The cost to pay is, as we shall see, less regular oscillations. 

A lower sensitivity to the very parameter choice is also seen in a generic model of 
pattern formation [16] , where it has been demonstrated that so called "quasi-patterns" , 
very similar to Turing patterns p2], occur even without fine-tuning of the parameters. 
Turing patterns are observed in a variety of systems, so along with them also quasi- 
patterns may be observed in very different applications, ranging from ecological systems 
as studied in [3j [16] to populations of players, playing a kind of prisoner's dilemma game 
[18] to neural networks [19] to genetic networks (7J [8j [9] . 

Similar observations have been made for systems more closely related to our unit. 
For example, quasi-cycles in the stochastic brusselator model [20] were detected well 
away from the regime of oscillations in the deterministic model. However, in contrast 
to these recent findings which could be well understood by treating demographic noise 
as a perturbation (although often amplified by stochastic resonance) on top of the 
deterministic behaviour, we shall see that noise has more dramatic effects on our system, 
going beyond the perturbative regime. 

The paper is organized as follows. In section [2] we present the deterministic and 
the stochastic formulations of the model. Section [3] deals with results of Gillespie 
simulations over different parameter ranges and as function of the system size. We 
show that quasi-cycles can be distinguished from limit cycles via a distinct, more rapid 
decay of their autocorrelation function, but apart from that they share other common 
features, particularly in the transition region. In section H] we shall solve the master 
equation approximately by using the van Kampen expansion in 1/Nq, where N stands 
for the system size. We derive the variances, autocorrelation functions and the power 
spectrum for the fluctuations of the two species, both in the fixed-point phase and in 
the limit-cycle phase. We show that the van Kampen expansion predicts the behaviour 
of all these quantities quite far from the transition points, where large excursions in 
phase space, induced by large fluctuations, are so rare that they can be neglected, but it 
breaks down close to fixed-point/limit-cycle transitions with large fluctuations and large 
excursions in phase space, which are characteristic for the excitatory behaviour of our 
model. This shows that finite-size effects are even more important in our model than 
in previously considered models of gene expression, and taking only first-order terms of 
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the expansion in 1/Nq largely neglects a characteristic feature of the BFU, that is its 
excitability. Conclusions and an outlook are presented in section |5j 




Figure 1. Motif of a bistable frustrated unit (BFU). Pointed arrows denote activation 
(increase in the production rate), the blunt arrow denotes repression (decrease in the 
production rate). 



2. The model 

Let us consider a simple model of frustrated bistability as depicted in Fig. [TJ In this 
model, A and B are two different species. Although from a statistical physics point 
of view the identity of A, B is not important, we shall assume here that A, B are 
two different protein types. This is in line with the original motivation of this model 
as a coarse-grained description of some genetic circuits. The protein A activates its 
own production (transcription in the biological language) and also the production of 
the protein B, which in turn represses the production of A. In this way, we have a 
self-activating bistable unit coupled to a negative feedback loop. As stated in the 
introduction, we shall call the motif from Fig. [TJa bistable frustrated unit (BFU). This 
motif exists as a part of many genetic circuits, such as the embryonic cell-cycle oscillator 
[2Tj and the circadian clock [15]. The simplest, idealised implementation of the BFU has 
been analysed on a deterministic level, with protein concentrations as the only dynamical 
variables, and it is known to produce oscillations in a certain range of model parameters 
[22] [TO] . In particular, the model studied in Ref . [10] , which we will adopt as our starting 
point, assumes that the protein production rates depend on the concentrations 0a and 
4>b of the two protein species as follows: 



a 



dt l + <f) B /Kl + 



(1) 



^f = 7(0.-0B), (2) 

where 7 is the ratio of the half-life of A to that of B. Here we shall focus on the case 
7 <C 1, that is when the protein B has a much longer half-life than A with a slow 
reaction on changes in A, while A has a fast response to changes in B. The parameter 
7 plays a similar role to the e-parameter in a FitzHugh-Nagumo element that separates 
the time scales of a fast variable, corresponding to the voltage variable in the context 
of neurons, and a slow variable, corresponding to the recovery variable there. The 
parameter K sets the strength of repression of A by B. We shall assume K 1, so 
that already a small concentration of B will inhibit the production of A. The parameter 
b determines the basal expression level of A. We set this parameter to anything larger 
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than zero but much smaller than one, so that the system cannot be absorbed in the 
state <pA = 4>b — and, simultaneously, the production rate of A is small for <pA ~ 0. 
In units where the production rate of B is equal to its degradation rate, K plays the 
role of a Michaelis constant that sets the strength of the repression of A by K. As 
argued in [10], the choice of the Hill coefficient (here as h = 1 in (<j)B/K) h ), used in the 
repression of A by B, is not essential. For the other Hill coefficients (activation of B 
by A with Hill coefficient h = and self activation of A involving powers of h = 2) we 
made the same choice as in [THl EE] for better comparability. The parameter a is the 
maximal rate of production of A for full activation (<j) A ^> b) and no repression ((f) b ~ 0). 
This will be our tunable parameter which we will use to control the behaviour of our 
model. This parameter seems to be the easiest one to control in real, experimental 
systems [10]. In the deterministic description of [13] . we extended the parameter range 
of a independently of the biological relevance. As a function of increasing a we found 
for small a a fixed-point regime with excitable behavior, followed first by a limit-cycle 
regime with an unstable fixed point and rich oscillatory behavior in this intermediate 
range of a, and next by a second fixed-point regime with excitable behavior for large 
a. The different regimes (or phases) are separated by subcritical Hopf bifurcations (for 
a definition see for example Ref. [23]) with corresponding hysteresis effects. Typical 
trajectories in phase space for the three regimes in the deterministic case are shown in 
Fig. 2 of [33] (as well as in Fig.s |2j [3] below), and the bifurcation diagram is displayed 
in Fig. 4 of the same reference. 

It is not at all obvious that Eqs. (D~H2]) can actually be interpreted as coarse-grained 
description of a full genetic circuit, in which other proteins [23] as well as mRNA [15] are 
typically involved. However, in the forthcoming paper [25] we shall show that Eqs. (CQ- 
[2]) can be derived from the assumed underlying reactions between genes, mRNA and 
proteins in a more detailed model, for a certain range of reaction rate constants. 

We now formulate an effective, stochastic counterpart of the model ([TH2J) by 
introducing individual molecules of A and B. At each time, the system is characterized 
by the number of molecules of each species, Na and Nb- The concentrations 4>a,4>b 
become discrete numbers: 4>a = Na/N and (fis = N B /N , where N plays the role of 
the system size. The parameter Nq is not the physical volume, but it allows us to control 
the average numbers of molecules A and B, which are proportional to Nq, and hence 
to investigate the role of demographic noise in the model. Molecules of both protein 
species are created and annihilated with certain rates. Since we want to establish the 
correspondence between the deterministic ([TH2]) and the stochastic version of the model, 
we assume the following rates for four possible processes: 



production of A 
decay of A 
production of B 
decay of B 



N A -> N A + 1 
N A -> N A - 1 
N B -> N B + 1 
N B -> N B - 1 



with rate N f(N A /N , N B /N ) (3) 

with rate N A (4) 

with rate jNa (5) 

with rate jNb (6) 
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where "rate" denotes the transition rate of a specific process, and f{4>A, 4>b) is given by 



It is by no means the only possible way of modelling the BFU on the level of individual 
molecules and chemical reactions, but we shall see later that this choice of the rates 
reproduces the deterministic system in the limit Nq — > oo. Also, this is only an effective 
model which neglects gene states, mRNA concentrations etc., which should be included 
in a more detailed model [25]. Similarly to the deterministic case (IlH2]), here we focus 
only on the dynamics of the proteins A, B. It should be noticed that the difference of 
half-lifes of species A and B here is still implemented via the parameter 7 in Eqs. (jSJ) 
and ((6]), making the effective production and decay rates of B much slower than those 
of A, without resolving the underlying genetic and mRNA-levels of how these rates are 
generated. 

The time evolution of our system can be described by the following master equation 
for the probability P(N A , Nb) for finding proteins of type A and Nb proteins of type 
B at time t: 



On the right- hand- side we have as loss terms for P(N A , Nb) the production of A with 
rate A /(Aa/A , Nb/N ), the decay of A proportional to Na, the production of B 
proportional to ■jNa, and its deletion proportional to jN B - Gain terms to P(Na,N b ) 
result from the decay of A out of a state with Na + 1 proteins, its production from a 
state with Na — 1 proteins with rate Nof((NA — 1)/Nq,Nb/Nq), the decay of B with 
rate ^(Nb + 1) and its production with rate jNa from Nb — 1 proteins of type B. 

The master equation (JHJ) is our starting point. In the limit of No — > 
00, our stochastic description becomes equivalent to the deterministic model 
E]) • One can see this by multiplying both sides of Eq. (jSJ) by Na or Nb and 
summing over Na and Nb- If we now define the average of some observable 



0(N A ,N B ) as (0(N A ,N B )) = En a ,n b °( N a, N b )P(N a , N b , t) and assume that 
(f(N A /N , Nb/N )) = f((N A )/N , (Nb)/N ) (which is fulfilled for a sufficiently sharply 



peaked probability distribution), we obtain the following equations for the averages 



f((j)A,(pB) 



a b + 4> 2 A 



(7) 



dP(N A ,N B ) 
dt 



= - (N f(N A /N , Nb/Nq) + N a + jN a + jN B )P(N A , N B ) 

+ (N a + 1)P(N a + 1,Nb) 

+ N f((N A - 1)/N ,N B /N )P(N A - 1, Nb) 

+ 7 (A B + l)P(N A , N B + 1) + lN A P(N A , N b - 1). (8) 



(N a ),(Nb): 




(9) 



d(N B ) 



j({N A ) - (N B )). 



(10) 



dt 
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Substituting (Na) = N <f>A, (Nb) = No4>b we arrive at our deterministic set of equations 
([I]-[2]) for the concentrations <Pa,<Pb- 

3. Computer simulations 

The master equation (J5J) cannot be solved exactly due to the strongly non-linear term 
f a-, <t> b) ■ Before we proceed to its solution by approximate methods, we shall discuss 
the results of computer simulations. We have simulated our stochastic model using the 
Gillespie algorithm [26], which samples the exact probability distribution P(Na, Ns,t) 
by simulating trajectories (iV^i), Ns(t)) with transition rates ()3H6]). The algorithm 
actually generates triples of points {Na,i, Nb,%, U}, where ti is the physical time at 
i-ih step of the algorithrqj- In this way, we obtain not only P(Na, N B ,t) for any 
finite time t, but also autocorrelation functions and power spectra of A^i(i), Ns(t) 
from the corresponding time series. For example, to obtain P(Na, Ns,t) from the 
Gillespie simulations, we determine a histogram that counts how often the trajectory 
visits a point (Na, Nb) at a given fixed physical time t. The stationary distribution 
P*(Na, Nb) = P(Na, Ns,t — > oo) is obtained by simulating the system for some time 
to until it reaches the stationary state, and collecting the histogram of visited points for 
all t > to- Probability distributions, derived from the Gillespie trajectories in this way, 
agree very well with direct numerical integration of the master equation (JSJ). 

In Fig. |2] we show deterministic (black lines) and stochastic (red symbols) 
trajectories obtained by numerically solving the deterministic equations ([IH2]) and 
via Gillespie simulations, respectively. Similarly to our previous work |13j . we use 
a as the control parameter, with other parameters kept fixed. Typically, we set 
K = 0.02,6 = 0.01,7 = 0.01 if not stated otherwise. As we vary a, we go from the 
deterministic fixed-point regime for a < ct\ ~ 31.10 (Fig. EK, d), through the limit-cycle 
regime (Fig. [2b, e), to another fixed-point regime for a > a 2 ~ 98.93 (Fig. [2b, f). In 
these simulations, after some time the system seems to converge to a neighbourhood 
of the deterministic fixed point and the deterministic limit cycle, respectively. The 
stochastic trajectories roughly follow the deterministic ones, with larger fluctuations for 
smaller system sizes (as seen from Figj2k-c, as compared to FigfJH-f). Although it is an 
abuse to speak of a fixed point and a limit-cycle trajectory in the stochastic system, we 
shall use these phrases to denote regimes of a which exhibit fixed-point or limit-cycle 
behaviour for N — > oo. 

If we, however, wait long enough, we see cycles also deeply in both fixed-point 
regimes, a = 15 and a = 150, see Fig. [3J top. As previously mentioned, we shall call 
such cycles "quasi-cycles" due to their absence in the limit No — > oo. It should be 
noticed that a in Figj3] bottom is close to, but still below the value of the deterministic 
bifurcation point and accordingly the deterministic trajectory converges to a fixed point, 

| The physical time t and the computer time i measured in Gillespie steps are not the same, nor are 
they strictly proportional, because natural time intervals between succeeding reactions are smaller or 
larger depending on the total rate of all reactions which changes with Na, Nb- 



Stochastic Description of a Bistable Frustrated Unit 



8 




300 600 900 1200 2000 4000 6000 4000 8000 12000 



(d) N A (e) N A ( f ) N A 

Figure 2. Trajectories in phase space up to Tq = 3000 (a) , Tq = 10 5 (b-c) and 
Tq = 5 x 10 5 (d-f) steps of the Gillespie algorithm, respectively, for Nq = 100 (a-c) 
and Nq = 1000 (d-f), in the fixed-point phases a = 15 (a, d), and a = 150 (c, f), and 
in the limit-cycle phase a = 50 (b, e). The initial condition is iVyt(i = 0) = 1000, 
N B (t = 0) = 100. 




50 100 150 200 250 300 350 400 450 500 



N A 

Figure 3. Trajectories in phase space for Tq = 10 4 ,a = 15 (top left), and 
Tq = 4 • 10 5 ,a = 150 (top right). Bottom: for a = 26, close to the transition region, 
we see cycles for the stochastic simulations (red points) and a trajectory (blue line) 
converging to the fixed point as deterministic solution. A*o = 100 in all cases. 
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but the stochastic trajectories show large fluctuations and proceed along cycles. The 
large fluctuations are typical for being in the vicinity of a transition point to another 
phase. The figure illustrates that it is not possible to visually distinguish the origin 
of the cycles on the basis of these trajectories. Similarly to finite-size shifts of the 
critical temperature in phase transitions, one interpretation here would be that the 
bifurcation point from the fixed-point to the limit-cycle regime is shifted for finite No, 
so that the observed cycles are genuine limit cycles in the limit-cycle regime at finite 
iVo- Alternatively, the stochastic system may be in the fixed-point regime, but with a 
large number of quasi-cycles. Therefore we need alternative measures to disentangle the 
origin of these cycles that we shall consider below. 

In Fig. HI left, we plot the frequency of such cycles determined in simulations as a 
function of a, for different system sizes Nq. The experimentally determined frequency 
agrees with the frequency of oscillations of the deterministic model in the limit-cycle 
regime, but it is clearly non-zero for a-values which belong to fixed-point regimes^. The 
smaller Nq is, the higher is the frequency of quasi-cycles in these regimes. For fixed a, 
the frequency decays exponentially with increasing Nq, which can be seen in Fig. 01 right, 
where we have plotted the frequency for a = 24 as a function of No. This suggests that 
quasi-cycles are large excursions from the fixed point, caused by crossing the inherent 
"energy barrier" which is characteristic for the excitable unit. 

In Fig. Owe show the deterministic flow ((INa/oH, dNs/dt) along with an example 
of a trajectory which ends in the fixed point. It shall illustrate that the height of the 
barrier to escape the fixed point increases with decreasing 7. The plot (left figure) 
for small 7 shows that (in almost all directions, apart from that towards decreasing 
Nb), the system can escape the vicinity of the fixed point only if a sufficiently large 
fluctuation drives it to the regime of small iV^/large Na, where the vectors of the flow 
point outwards from the fixed point. For an excitation beyond such a large barrier the 
system makes a large excursion in phase space, while it directly relaxes to the fixed 
point if the excitation is below this threshold. On the other hand, small N values help 
to cross the barrier, because the fluctuations are then larger. The plot for larger 7 (right 
figure) illustrates that in this case already small fluctuations are sufficient to escape the 
fixed point, such fluctuations will then lead to small excursions (cycles) in phase space. 

We have also measured the number of quasi-cycles as a function of 7, other 
parameters being fixed (results not shown). The increase of the height of the barrier 
with decreasing 7 manifests itself in an increasing number of quasi-cycles when 7 grows. 

We shall show that these cycles, interpreted as quasi-cycles, can occur as quite 
regular oscillations close to the transition region. In Fig. (6J top left, we plot NA{t) 
for a = 28 and Nq = 10 3 . In the same figure, top right, we plot the power spectrum 
Ppf A (oo) = J Nji(t)e~ ltuJ dt\ which shows a peak at u = 0.027. This clearly indicates 
the regular, oscillatory nature of quasi-cycles for a close to the transition point. As a 
gets smaller, these oscillations become less regular and finally become rare stochastic 

§ A similar effect has been observed in Ref. 1 1 3 j where multiplicative noise is introduced in the 
deterministic equations ([T][2]). 
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a N 



Figure 4. Left: Frequency of (quasi-)cycles as a function of a. Red circles: iVo = 10 3 , 
blue squares: iVo = 10 4 . Black line: frequency of oscillations (inverse of the oscillatory 
period) of the deterministic model. Right: frequency for a = 24 and different Nq. 
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7 = 0.5 



10 



15 



20 



Figure 5. Plots of the deterministic flow (dNA/dt,dNB/dt) = No(d(f>A/dt,d(/>B/dt) 
from Eqs. (Tj][2]) (blue vectors) together with an example of a deterministic solution (red 
curve) in the fixed-point regime a = 15. For 7 = 0.5 small fluctuations in all directions 
from the fixed point (end point of the trajectory) lead back to the fixed point after 
small excursions (oscillations) in phase space, while for 7 = 0.01 small fluctuations in 
almost all directions would directly lead back to the fixed point without any excursion, 
so that the barrier for an excursion is very high; only for a small fluctuation in the 
direction of negative Nb a large excursion would be possible. Correspondingly, the 
vector field is almost constant over large areas, but changes very rapidly over small 
areas, in contrast to the flow for 7 = 0.5. 



events, as shown in the same figure, bottom. 

In the following we shall show how we can distinguish genuine limit cycles from 
quasi-cycles by comparing decay rates of autocorrelation functions of N A (t), iYs(t) 
obtained from simulations. For simplicity, we shall focus only on iV^. In Fig. [7] we 
plot the autocorrelation function: 

c ( ) = fo™(N A (t) - (N A ))(N A (t + r) - (N A ))dt 

J tin -(N A (t)-(N A }ydt 

Here t max denotes the length of the sample in time (optimally, as long as possible), and 
the value of N A (t) obtained from the Gillespie simulation is such that N A (t) = N A (ti) 
for the largest ti such that tj < t. We stress that C A {r) measures autocorrelations in 



Stochastic Description of a Bistable Frustrated Unit 



11 



4000 
3000 

<; 2000 



1000 








fe; 2000 



-JJUI J Ji 



1000 2000 3000 4000 5000 




1000 2000 3000 4000 

t 



5000 




Figure 6. Plots of N^it) (left) and their power spectra (right), for iVo = 10 3 and 
a = 28 (top) and a = 20 (bottom) in the fixed-point regime. Quasi-regular oscillations 
are visible for a = 28. For a = 20, the system still shows spikes in Na, but they are 
much less regular. 



< 
O 




100 200 300 400 500 600 

T 



< 

o 



1 




0.8 






a 


0.6 


H 







0.4 


\ 




■ 


0.2 






\ 









-0.2 





a=150 
N =100 



200 400 600 800 1000 



< 
O 



1 

0.8 
0.6 
0.4 
0.2 


-0.2 



1.2 

\ a=5 ° 0.8 
| N =100 Q 4 








WW 






i o 








\ f% 250 500 7501000' 


I I 


1 


" mi if » 
1 



200 400 600 800 1000 



Figure 7. Autocorrelation functions Ca(t) for A^(t) for different values of a = 
20, 50, 150 and No = 100. The inset for a = 50 shows the deterministic behaviour 
of Ca{t) in the limit-cycle regime. The maxima for r > indicate the presence of 
(quasi-)cycles. 
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Figure 8. Decay rate /3 of the autocorrelation function of Na as function of a, for 
N = 100,7 = 0.01. The rate /3(a) changes its slope at a = 26 ± 6 and a = 96 ± 9 
which agrees with the lower and the upper bifurcation point ai ~ 31.10, ot2 ~ 98.93. 



the real, physical time and not in the number of steps of the Gillespie algorithm. In 
practice, we have calculated Ca{t) from the numerically obtained time series {iV^ij), ij} 
by resampling them at uniform time intervals, performing the Fast Fourier Transform 
(FFT), squaring the absolute values of the coefficients, and doing the inverse FFT. This 
method is not only much faster than calculating Ca(t) directly from Eq. ( iTTTl . but it 
also produces the power spectrum P{ui) of the signal as a by-product, which we also use 
later in this work. Figure [7] shows examples of Ca{t) for all three regimes of a. One sees 
a number of maxima for r > which indicate the presence of cycles and quasi-cycles. 
The amplitudes of these peaks decay exponentially with r. If we fit the exponential 
decay e _/3r to these peaks and plot the decay rate (3 as function of a as in Fig. El we 
see that (3 changes faster with a for quasi-cycles than for limit cycles. This behaviour 
seems to be a generic feature of oscillating stochastic systems [161 EZ]- If we now fit 
straight lines to (3(a) in different regimes (see Fig. [8]), we see that these lines cross very 
close to the deterministic transition points 01,02. 

We shall now show that these transition points, which were identified as subcritical 
Hopf bifurcations in the deterministic limit, have precursors in the stochastic model. 
In general, such precursors are larger fluctuations, coherent behaviour in space and/or 
time, and phenomena similar to critical slowing down known from second order phase 
transitions. In Fig. [9l we plot sets of 10 independent trajectories of N a{£) for iVo = 10 5 , 
7 = 0.5, starting from the same initial state, each set for another value of a. We observe 
an increased coherence in time: as a approaches the transition at 0:1(7 = 0-5) ~ 45 from 
the fixed-point to the limit-cycle regime, independent time series follow each other for an 
increasing period of time. For example, the time series for a = 44 follow the very same 
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Figure 9. Plots of N A {t) for No = 10 5 ,7 = 0.5 and different a from the fixed-point 
regime: 25 (black), 40 (red) 43 (green) and 44 (blue), curves from bottom to top. For 
each a, 10 independent realisations of N A (t) starting from the same initial condition 
are shown. As a approaches the transition point a ~ 45 to the limit-cycle regime, 
the coherence time between independent realisations increases. This indicates that 
relaxation to the stationary state slows down at the transition point. 

oscillating trajectory up to T > 100, while away from the transition region (a = 25), 
trajectories reach the stationary state already after T 20. To observe the coherence, 
trajectories have to stay close to the deterministic solution. This is fulfilled if the number 
of particles is sufficiently large. For 7 = 0.5, it is sufficient to choose N = 10 5 . For our 
usual parameter choice of 7 = 0.01, it is necessary to take N > 10 6 to see the coherence, 
and even then the increase in the coherence time is seen only in a very narrow region 
below the critical value a\ ~ 31.1 (results not shown). 

In the last part of this section, we shall discuss the stationary probability 
distribution P*(N A , N B ) = P(N A , N B , t — >■ 00) of finding the particle at (N A , N B ) after 
the system has forgotten the initial condition. The plot is shown in Fig. [10] for two values 
of a, in the fixed-point and in the limit-cycle regime. The support of the probability 
distribution has always a donut-like shape and a maximum for small values of N A . The 
relative amplitudes of the donut and the peak depend on a. The peak corresponds to 
the probability localized around the fixed point, whereas the donut shape corresponds to 
cycles and quasi-cycles. We also see from Fig. [10] that the range of fluctuations strongly 
depends on the location along the limit cycle: it is very large for large values of N A and 
small for small values of N A . In the next section we shall analytically predict how the 
variance of N A evolves in time, and how the autocorrelations decay in time for different 
regimes of a. 

4. Analytical approach 

So far, we have shown the results of computer simulations. We cannot find analytical 
solutions of the master equation in a closed form, nor is a product ansatz P(N A , N B ) = 
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Figure 10. Density plots of the probability P*(Na, Nb) in the stationary state, for 
a = 15 (fixed-point regime, left) and a = 50 (limit cycle, right), for Nq = 100. Insets 
show peaks in the probability from small- Na regions. 



P(Na)P{Nb) justified, since both variables are not on an equal footing with respect 
to their coupling in our model. However, we shall see that the van Kampen expansion 
[28] in the inverse of the system size Nq captures many (but not all) of the observed 
features. This method assumes that we can split the integer variables Na, Nb into a 
deterministic part proportional to N , and a stochastic part proportional to y/No: 

N A (t) = NoMt) + V^Vof (*) (12) 

N B (t) = N <p B (t) + \^N~oV{t), (13) 

where fluctuations £ of Na and rj of Nb are continuous variables of order 0(1). We now 
introduce the step operators 

E+[g(n)} = g(n + 1), E-[g(n)} = g(n - 1), (14) 

acting on an arbitrary function g(n). The master equation (jBJ) can then be rewritten as 
follows: 

8P{Na,Nb) = (E _^ _ 1)NQf{NA/N ^ NB/No)P{N ^ NB) 

+ (E+ A + 7 E Wfl - 7 - 1)N A P(N A , N b ) 

+ (E+ fl -l) 7 iV B P(iV A ,iV B ), (15) 

where f(NA/N ,N B /N ) is defined in Eq. (J7J). For given, yet undetermined functions 
4>A(t), 4>B(t), Eqs. ( fT2TfT3l) correspond to a transformation of variables (Na,Nb) to 
(£,77). Therefore the function P(Na, N B ,t) defines a new function H(£,r),t). This way 
the time derivative of P(Na, Nb, t) can be written in terms of the new variables as 

dP dU r—d<P A dU r—d<f)BdU 

■W = ~dt " ^°lTf ~ * °~dT~dri' (16) 
If we make a Taylor expansion of the step operators in powers of Nq 1 ^ 2 , 

E ^ = l±^| + ^o- 1 ^ + - (17) 
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-i d 1 
1 ± iV 2 — + -N, 



drj 



drj'' 



+ ... 



(18) 



_____ ^ ^2 

and insert these expressions into the master equation (|15p . we obtain to order N Q the 



deterministic equations ([IH2]) and to order N® 

an(x,t) N - 

at " ^ ij dx{ 



1 the following Fokker-Planck equation: 



_^ an(x,t) 



(19) 



where x = (xi, x-z. 



and 



.4 



B 



i,ie{i,2} 

(£,77), and Ajj and .By are the following matrices 
+ /b 



7 







-7 





7(0A + <\> 



B 



df(4> A (t),<l>B(t)) 



(20) 



(21) 



and 



Here we have introduced the notation / = </>#(£)), /a — d<j> A {t) 

fs = d ^^Q^1ft)^ • This Fokker-Planck equation is known to have a Gaussian 
distribution as unique solution [2E], so it is sufficient to determine its first and second 
moments. From Eq. (fT9l we derive the following equations for the first moments 

d(0 



alt 
alt 



(Ja- 

7«e> 



i)(0 



f B {v), 



(22) 
(23) 



with (£) = J £n.(£,r),t)d£dr) and (77) = J r}H(£,r},i)d£dr}. For the second moments we 
obtain 



d 
~dt 







2(f A - 1) 





2/b 


(v 2 ) 







-2 7 


2 7 






7 


/b 


(fA-l 






<pA + f 








+ 


7(0A + 0b) 


















7) 



(24) 



We shall show now that these equations correctly describe fluctuations not only in the 
fixed-point regime, but also in the limit-cycle regime, provided that N is sufficiently 
large and we are not too close to the transition (bifurcation) points. We shall also show 
how this approximation breaks down in the vicinity of transition points, and explain the 
reason why the van Kampen expansion does not work in these cases. Therefore, we will 
investigate how the variances of Na, Ng evolve in time, calculate the autocorrelation 
function and the power spectrum in the stationary state in both regimes, and compare 
these quantities with Gillespie simulations. 
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Figure 11. Stationary-state variance (£, 2 ) s as a function of a, obtained from Eq. (f25|) 
(solid line) and from Gillespie simulations for N = 100 (blue circles), 10 3 (green 
squares), 10 4 (black diamonds) and 10 5 (red triangles). 



4-1. Solutions in the fixed-point phase 

Let us first discuss the solution for a from the fixed-point phase. Assuming that the 
system starts always from the same point in the phase space so that the averages 
(£(0)) = and (77(0)) = are zero at time t — 0, the linearity of Eqs. (|22H23|) causes 
the averages to stay equal to zero for any time t > 0. Therefore, the equations for the 
second moments are actually the equations for the variances. In the stationary state, 
the time derivatives on the left-hand side of Eq. (I24j) vanish, and <pA = <Pb = <ft*, where 
<fi* can be determined numerically from Eq. ([I]) as the root of <p* = f (<f)* , <p*) . We then 
obtain from Eq. fl24"j) that the stationary second moments of the fluctuations are 



where 



{e)g = + + (25) 

(v 2 )s=<P* + ^ 1 + il D fA)fB , (26) 
(ft). = 20* 7 + (1 ~ /A)/B , (27) 



D = 2( 7 + 1-/ A )(l-/ A -/ B ). (2* 



The subscript "s" stands for the stationary state. In this set of equations, the quantities 
4>Ai 4>Bi Ia = Sr an<1 fs = Sr should be evaluated at the fixed point which depends 



d<f> A 

on the choice of parameters , -f,b,K,a. In Fig. [11] we compare the variance (|25|) with 
(£ 2 )s = {(N\) — (Na) 2 )/Nq obtained from the Gillespie simulations. The agreement is 
very good if a is deep in the fixed-point regime, but it becomes worse as a approaches any 
of the transition points. This is caused by the more and more frequent occurrence of large 
quasi-cycles, induced by large fluctuations, not captured by the above approximation, 
which increase the variance by orders of magnitude. 

In a second step we would like to derive the autocorrelation function for £(t) from 
the solutions of the moment equations. The differential equations for the time evolution 
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Figure 12. Comparison between autocorrelation function Ca(t) from Eq. (thin 
lines) and the results of Gillespie simulations (thick lines). Left: No = 5000, a = 
15,20,25 (black/grey, red/orange, blue/turquoise). Right: N = 10 5 , a = 25,27,29 
(black/grey, red/orange, blue/turquoise). For the same value of a = 25, the agreement 
gets better when Nq increases. 



of autocorrelation functions (£(0)£(r)) and (£(0)7?(t)) are the same as for the averages 
( I22H23]) . since £(0) is constant from a r's point of view, but the initial condition must 
take into account fluctuations around the fixed point, which do not average out to zero 
when combined with £(r). We therefore have: 



air 



(/,-l)(^(r)) + / B PMr)), (29) 
7«£(0)£(t)> - <£(0Mt)», (30) 



dr 

with (£(0)£(0)) = (£ 2 ) s , (^(0)77(0)) = (£77), from Eqs. ([25]), ([27]) as the initial condition. 
In the following we focus only on the autocorrelations of £(t). As r-dependent solutions 
we obtain 

mar)) = ^ (e^ + e^) 

+ 2(Ax-A 2 ) [e e h {61} 



where A12 are given by 



1 - f/ 



- ± \ V(7 - !a + 1) 2 + 47(^ + /b-1). (32) 



2 2 

Again, 0^, ^-dependent quantities should be evaluated at the fixed point. In Fig. [T2] 
we compare the numerically obtained Ca[j) from Eq. (fTTj) with 

r , r ^ = (jww) - (iv^o)) 2 = mm m , 
A[ ) (Nim - (n a (o)) 2 (em ' [66) 

where (£(0)£(t)) is calculated from Eq. (|3T|) . Similarly to what we have seen for the 
variance, the range of a over which the above formula applies increases with Nq. 

The power spectrum of the fluctuations can be either obtained via the Fourier 
transform of the corresponding autocorrelations, or by directly Fourier-transforming 
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the Langevin equation for £(£),r/(£) that corresponds to the Fokker-Planck equation 
(fT9j) [28]. In our case, the Langevin equation is given by 

^£>=A*x(0 + CW, (34) 

where x = {C,,rj}, the drift matrix A* from Eq. fl20l) is evaluated at the fixed point 
(pA = 4>b = 0*, and C(t) = {(i(t), C2(t)} is a bivariate white Gaussian noise with 



at)) =0, (35) 

%m(t'))= Btjit-t), (36) 

where i,j G {1,2} and the matrix B* is given by Eq. (121j) evaluated at the fixed point. 
As before, we shall focus on £(£). We transform Eq. fl34l) to the Fourier space using 
fj(o;) = f 00 Xi(t)e~ luJt dt, solve the resulting algebraic equations for £,fj, and calculate 
the power spectrum -P^(w) = (\C,(u)\ 2 ), where the average can be expressed in terms of 
the coefficients B*j. Remembering that <Pa = 4>b = <P* and /(</>*,</>*) = <ft*, we finally 
obtain 

p (u) = 20- [ 7/ l + 7 2 + 

^ J ^ + [(l- /A ) 2 + 7 2 + 27/B]u; 2 + 7 2 (l-/A-/i 3 ) 2 ' 1 J 

where /a, Jb should be evaluated at the fixed point. A comparison of the results 
obtained from Gillespie simulations with the above result ( 157)) is displayed in Fig. [I5J 
Interestingly, the analytic formula predicts a peak in P^(oj) for a — > ax, which agrees 
with simulation results for iV large enough. This peak indicates the presence of quasi- 
cycles. However, when a approaches the critical point, the position and the shape 
of the peak obtained in simulations differ from that of Eq. (1371) . Quasi-cycles which 
correspond to a large excursion in phase space are triggered by large fluctuations. Such 
fluctuations are very rare deep in the fixed-point regime, and due to their non-Gaussian 
nature, the induced quasi-cycles are not covered at all by the van Kampen expansion. 
Close to the transition region, large fluctuations are naturally more frequent, along 
with them large excursions in phase space, so the discrepancy between the van Kampen 
results and the Gillespie simulations are more pronounced. On the other hand, the van 
Kampen expansion is able to predict the occurrence of quasi-cycles which are induced 
by small (Gaussian) fluctuations, leading to small excursions in phase space similarly to 
the quasi-cycles in the brusselator system [20J. Excursions of this type are more frequent 
for a-values away from the transition region; it is these quasi-cycles which induce a peak 
in the power spectrum at non-zero frequency both in the van Kampen expansion and 
the Gillespie simulations with a reasonable agreement. 



4-2. Solutions in the limit-cycle phase 

It turns out that the van Kampen size expansion works not only in the fixed-point 
regime, but also (to some extent as we shall see) in the limit-cycle regime, where 
the expansion (fT2] - [T3l has to be carried out about the time- dependent trajectory. 
We begin with calculating the variances (£, 2 {t)), (f] 2 {t)) along the limit cycle. Our 
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Figure 13. Comparison between the power spectrum P^(u) in the fixed point from 
Eq. (pT7]) and obtained in simulations for N n = 5000 (top) and N — 10 5 (bottom). 
Top left to right: a = 15, 20, 25. Bottom left to right: a = 27, 28, 29. 

starting point is again Eq. (jMJ), but now with time- dependent functions (pA(t), </>b(£), 
fA(t), /b(£), fit), which should be determined from the limit-cycle solutions of the 
deterministic equations. Unlike in the fixed point, the integration cannot be done 
analytically in the limit-cycle regime. However, it is straightforward to integrate Eq. (T24|) 
together with Eqs. (JTH2D numerically, assuming some initial values of 0^(0), 0#(O), and 
(f 2 (0)) = (r? 2 (0)} = (£(0)77(0)) = 0. As before, the mean values (£(*)) = (rj(t)) = 
for any time t > 0. Altogether, we have five differential equations, which can be 
integrated by any method for solving lst-order ordinary differential equations. In 
Fig. [14] we compare the variance (£ 2 (t)) obtained in this way for a = 50 (deeply in the 
limit-cycle regime) with the variance calculated from Gillespie time series according to 
(£ 2 (t)) = ((A r y i(t) — N (l)A{t)) 2 )/N , where 4>A(t) is the deterministic solution (obtained 
numerically) and (...) denotes the average over many simulations starting from the 
same initial state 4>a{0) = O,0_b(O) = 0.5. The variance oscillates with the same period 
T(a = 50) ~ 178.1 as the deterministic solution 4>A{t), 4>B{t)- The variance changes by 
many orders of magnitude along the limit cycle, with sharp peaks corresponding to fast 

changes in <pA{t), see Fig. [TH right. For N = 10 5 , the quantity \J /4>A(t) which 
measures the ratio of fluctuations to deterministic concentrations becomes of order 1 
only for a brief period of time. For smaller No = 10 3 , however, this ratio is larger than 
one for a significant part of the oscillation period, which means that the van Kampen 
approximation loses its applicability. As expected, the agreement between theory and 




Figure 14. Comparison between the variance of £(t) as function of time, obtained from 
100 independent Gillespie simulations (grey line), and via the van Kampen expansion 
(black line), for a — 50 and Nq = 10 3 (left) and 10 5 (middle). The variance strongly 

varies along the limit cycle. Right plot: comparison of the fluctuations J ^ ^ (black 
line) relative to the concentration (f>A(t) (green thick line) for Nq = 10 5 . 



simulations is thus better for large Nq, exactly as for the fixed-point regime, and the 
time period over which the van Kampen expansion applies also increases with Nq. 

Now, in principle we would like to find the autocorrelation function and the 
spectrum of infinite time series (that is for t max — > oo), exactly as we did in the fixed- 
point regime, for finite Nq. This is the stationary limit of a stochastic system with 
non-negligible demographic fluctuations. In the appendix we shall show that (i) the van 
Kampen expansion is valid only for t max <C O(NqT), where T is the period of oscillations 
and so it breaks down for infinite time series, (ii) for very short t max (a few periods of 
oscillations), the deterministic limit iVo — > oo is actually a reasonable approximation 
also for Nq < oo and the van Kampen expansion is not really needed, and (iii) in the 
physically relevant limit (Nq < oo,t max — > oo) the results of the Gillespie simulations 
can be analytically reproduced if we take into account that the physical time t itself is 
a random variable, which is neglected by the van Kampen expansion. The randomness 
of t results from the discreteness of our system which changes its state in discrete time 
steps of various length. This is taken into account in the Gillespie algorithm where 
time is incremented at each step by an exponentially distributed random number. As 
explained in the appendix, we model this by assuming that variables such as iV^, jVg 
are periodic functions of some new time x that plays a similar role to the number of 
Gillespie (time) steps. For a given fixed physical time t, the number of Gillespie steps 
varies due to random reaction times. Similarly here, for a fixed time t also x varies 
due to stochastic increments W(t), which shall reflect the fluctuations in reaction times. 
More precisely, x is assumed to describe a Wiener process (Brownian motion) according 
to dx/dt = 1 + W(t) with random numbers W(t) following a Gaussian distribution 
with variance a 2 . It should be noticed that we assume here the increments W(t) to be 
the only source of stochasticity, while we neglect demographic fluctuations in iV^, Nb 
and also the possible dependence of time fluctuations on the position along the limit 
cycle. As explained in the Appendix, this is valid for "spiked" oscillations which exist 
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Figure 15. Spectrum of NA(t) for N = 4000,7 = 0.01, and 2 26 Gillespie steps (grey 
circles), averaged over 100 realisations, compared to Eq. ((591) (black line). 



in our model for 7 -C 1. We then calculate the power spectrum by averaging over these 
fluctuations in x: 



N A {UJ) 



dt 



dt'e^-^lM^Mxit')) 



(3* 



where the average is over different positions in the new time x(t), x(t') of the Brownian 
motion at physical times t,t', and 4>a{%) is the same deterministic solution as before, 
but in the time variable x. For more details we refer to the appendix. The result for 
the power spectrum is 



\f n \ 2 [n 2 a 2 (n 2 u 2 + n 4 a A /4 + u 2 )} 

(n 4 a 4 



^ 16n 4 c4 + 8n 2 u 2 (n 4 a 4 - 4co 2 ) + (n 4 a 4 + 4u 2 ) 2 ' ^ 

where f n are Fourier coefficients of </u(i). We neglect the zero mode ~ / 2 u" 2 because 
it can be removed by shifting the average N A to zero. Equation (15^]) shows that each 
Fourier mode n, which would give a Dirac-delta peak for N — > 00, is smeared out to 
a skewed Lorentzian-like function for iVo < 00. In Fig. rT5]we compare this result with 
exact Gillespie simulations for iVo = 4000,7 = 0-01 an d with a 2 fitted to the data. The 
agreement is impressive and it suggests that we have indeed captured the main source 
of randomness in Gillespie simulations: the fluctuations in reaction times rather than 
in N A , N B . 



5. Conclusions and Outlook 



We have considered the fully stochastic version of a bistable frustrated unit which 
consists of a self-activating loop that activates its own repression. The unit is a simplified 
model of many genetic circuits. In the deterministic limit this unit shows excitable or 
oscillatory behaviour depending on the parameters. In our stochastic model, we have 
found quasi-cycles also deeply in the regions which correspond to excitable behaviour 
in the deterministic limit. The smaller the system size, the larger the fluctuations, 
and the more frequently the quasi-cycles occur. Therefore, in comparison to the 
deterministic version of this unit one may conclude that less fine-tuning is needed to 
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sustain oscillations. On the other hand, quasi-cycles can still be distinguished from true 
cycles by looking at the autocorrelation function decay rate or the power spectrum. 

It is known for systems such as deterministic FitzHugh-Nagumo units that their 
behaviour is much more versatile when these units are coupled with delay in the excitable 
regime than in the oscillatory regime [29]. It would be interesting to see whether coupling 
of BFUs, of which we considered a deterministic version in Ref . [13] , produces similar 
effects and in particular, whether the quasi-cycles of our single unit remain dynamically 
equivalent to genuine limit cycles of the oscillatory regime when these units are coupled. 
This will shed some light on the question of robustness with respect to a suitable 
parameter choice to maintain oscillations. 

In natural genetic systems it may be challenging to disentangle oscillations which are 
due to internal fluctuations from those with genuine limit cycles. What is the "normal" 
mode of performance in these systems? Are quasi-cycles an adequate replacement of 
limit cycles even in case of delayed interactions? Most likely there is no universal answer 
to this question, but it will depend on the very system. 

A large variety of natural genetic circuits have more complicated bistable systems 
coupled to negative feedback loops. Examples for such systems are the cAMP signalling 
system in the slime mold Dictyosthelium Discoideum [30] , the embryonic division control 
system [3TJ[32], or the MAP K- cascade [33]. In all these systems the number of reacting 
participants as well as the reaction events fluctuate, so that one should be aware of 
possible effects whose right interpretation will be in terms of inherent fluctuations. 

Lastly, in Sec. 4.2 we have developed a new method of calculating power spectra of 
oscillating stochastic systems, simulated with the Gillespie algorithm and corresponding 
to biochemical reactions. The method shows that the most important source of 
stochasticity is the stochastic nature of the time variable t and not the demographic 
fluctuations £,r) about the classical trajectory of concentrations ipAi'pB- It would be 
very interesting to see if this remains true also for other systems like the stochastic 
brusselator model [20] . 
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Figure 16. (£ 2 (i)) obtained fr om the van Kampen expansion grows linearly in time. 
Results for a = 50. 



6. Appendix 

(i) Let us first explain why the van Kampen expansion breaks down for large times. We 
are interested in calculating the spectrum 



The above two-point correlation function has an O(Nq) contribution from the 
deterministic equation, which is periodic in t,t', and an O(Nq) contribution from 
stochastic fluctuations. These fluctuations have also periodic oscillations in t, t', they are 
however superimposed on a linear growth, see Fig. [16] Due to this unbounded growth, 
the term Nq eventually overcomes the term NQ<pA{t)4>A(t') for sufficiently large 

t, t' . This is physically unrealistic; real fluctuations cannot depart that much from the 
classical trajectory, because they are at most of order v% for large Nq unless the system 
is close to a transition or bifurcation point. As already noted in Ref. (20], this limits 
the applicability of the van Kampen expansion to short time series t max <C O(N T), 
and hence the method cannot accurately predict the spectrum P^ a (uj) for infinite time 
series as we want. 

(ii) Before we proceed to an alternative method of finding PN A (t) for t max — > oo, let 
us consider t max « a few periods T. As mentioned above, the expansion fll2|13j) works 
in this regime, but its contribution to the spectrum is ~ 1/Nq of the contribution from 
deterministic 4>a{P)- We may thus hope that the deterministic spectrum P$ A {t) should 
be already a good approximation to the exact, finite-size P^ A (u). This is indeed true 
for iVo large enough. To see this, let us calculate the Fourier transform of <pA{t) for a 
finite period of time t max : 




(40) 



Using Eq. ( !T2l we can express (Na^Na^')) as 

(N A (t)N A (t')) = N*Mt)M*) + (atw)) ■ 



(41) 




(42) 
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Figure 17. Left: plot of the spectrum P/v A (w) from Gillespie simulations (red points, 
averaged over 100 realisations for Nq = 10 5 ,7 = 0.01) and Eq. (|45|) for deterministic 
4> A (t) (line). Initial condition was N A (0) = 0.0354896^,^5 = 0.5190827V which lies 
on the deterministic limit cycle, and the time series of length t max = 1904 corresponded 
to 10.7 periods of oscillation or about 2 27 Gillespie steps. Right: comparison between 
numerically obtained spectra Pn a (oj) for iVo = 4000,7 = 0-01 an d different number 
of Gillespie steps (different lengths of time series) 2 22 , . . . , 2 26 (curves from purple to 
black) . 



For simplicity, let us assume that initial values 0a(O), 0b(O) lie on the limit cycle. Then, 
4>A{t) is periodic (with period T) and it can be written as a Fourier series, 



<t>A{t) = E f nt 



i2irnt/T 



(43) 



The coefficients /„ can be determined numerically from </>a(£)- Inserting Eq. (143]) into 
Eq. (l4"2l) and performing integration we obtain 

i E ~ o-_/^ . ( 44 ) 



0J 



2vm/T 



so that the spectrum is finally given by 



(45) 



The result will generally depend on the initial condition 0^(0). This must be taken 
into account when comparing with simulations which should be performed for the same 
initial condition as the analytical calculations: Na(0) = No4>a(0), Nb(0) = A/o</>b(0). 
In Fig. [TTl left, we compare Eq. ( l45l) truncated to the first 30 terms in the Fourier 
series with Ppj A (t) from Gillespie simulations for short t max ~ 10T. The agreement is 
quite impressive, especially when taking into account that we have completely neglected 
stochastic effects in our calculations. Amplitudes and widths of the peaks are predicted 
correctly, the only discrepancy is seen in valleys between the peaks, where stochastic 
terms give significant contributions. 

The above method cannot be applied to longer times, in particular we cannot hope 
that it will give correct results for t max — > oo. In this limit, stochastic effects become 
significant. This is seen in Fig. [T71 right, where we plotted the spectra from simulations 
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for an increasing length of time t max . We see that the amplitudes/widths of the peaks 
first increase/decrease with increasing t max , but soon they stabilise at some finite value. 
However, according to formulas f!44p and ( )45|) . the peaks should become Dirac-delta 
functions for t max — > oo. Obviously, the deterministic approach fails in this limit, and 
so does the van Kampen expansion. In the rest of this section we shall show how to 
avoid these problems and how to reproduce the spectrum measured in simulations. 

(iii) We have already shown that stochastic fluctuations of give only a small 
contribution to Pn a (uj) for large Nq. Now we shall argue that the main contribution 
to the spectrum, which is responsible for the broadening of peaks for any finite No, 
comes from stochastic fluctuations of the physical time t. Although this is neglected 
in the van Kampen expansion, in stochastic systems of discrete particles t is not a 
continuous variable but it advances in discrete steps, from one event (change in the 
variables Na, Nb) to another one. The length of these steps depends on the position of 
the system in phase space. This feature is captured in Gillespie simulations, in which t 
is a random variable, incremented with each step of the algorithm by an exponentially 
distributed random number. 

To model the stochasticity of t, we shall introduce another variable x which could 
be identified with a properly rescaled Gillespie time (number of steps) such that x(t) is 
a Wiener process (Brownian motion) with drift and it evolves according to 

dx/dt = 1 + W{t), (46) 

where W(t) is Gaussian noise with zero mean and variance cr 2 < 1 which controls the 
strength of fluctuations (and thus it should decrease with the increase of No). In the 
limit cr 2 — > we have t = x that corresponds to the deterministic limit. The variance 
cr 2 should in principle depend on x, but here we assume it to be constant and equal to 
the average value taken over one period of the limit cycle. This assumption simplifies 
calculations very much and, as it will turn out, is a very good approximation. The reason 
is the following. We are mainly concerned with the case 7 <C 1, for which oscillations in 
iVyi(t) are "spiked", i.e., they are narrow peaks rather than smooth, sine-like oscillations. 
To be definite, let us consider fluctuations in the physical time which has elapsed after 
the system has performed one period in phase space. The main contribution to the 
stochasticity in this time comes from the short duration of the peak in iV^, where the 
number Na changes drastically over a few Gillespie steps (a rapid increase followed by 
a rapid decrease), whereas many reactions take place in the region where Na changes 
slowly, so that fluctuations in the number of Gillespie steps there average out, the better, 
the more Gillespie steps, and the larger the system. This means that the fluctuations 
in the physical time for a full period are dominated by the temporal duration of the 
peak in NA(t) in relation to the number of Gillespie steps during this peak: the sharper 
the peak, the lower the number of Gillespie steps, the stronger the fluctuations of the 
physical time when the peak occurs. For large iVo the height of the peak, i.e., the 
maximum number of Na, is proportional to iVo and stays almost the same from cycle 
to cycle (corresponding to small demographic fluctuations for large Nq); so it is just 
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the phase of N^t) at which the peak occurs within the limit cycle which fluctuates 
from cycle to cycle. The resulting fluctuations in the physical time for a full period 
average out more slowly (since they are less frequent) as compared to the demographic 
fluctuations and the fluctuations in the Gillespie steps outside the peak. The latter 
two behave similarly and decrease with increasing system size. The overall number of 
Gillespie steps is proportional to N . The variance a 2 should therefore decay as ~ l/iV , 
as a result of summing up O(N ) exponentially distributed random variables, each of 
them with variance 0(1/Nq), similarly to the demographic fluctuations, so in principle 
both effects would give the same contribution of 0(l/y/No) to A^(t). This allows us 
to neglect the contribution from demographic noise and noise in Gillespie steps outside 
the peak, and to concentrate on the stochasticity in time due to the peaks when they 
dominate the fluctuations. 

We shall now assume that the evolution of Na = Nq(()a is fully deterministic in this 
new variable x and that the only randomness in the evolution of N A {t) is due to the 
stochastic nature of t(x). The spectrum P Na (u>) is then given by 

P Na (co) = N 2 r dt f°° d^-^U A {x{t))<f> A {x{e))) . (47) 
Jo Jo N ' x 

The average ^ . . . \ is over different realizations of the Wiener process x(t) and 4>a (x) is 
the same deterministic solution as before, but in the new variable x. Using the Fourier 
series representation of </>a{%), we have 



oo 



M*(t))M*(f))) = E E fnfm{e iiJ0[nxit)+mxit ' )] ) , (48) 



n=— oo m=— oo 



where we have introduced uq = 2ti/T. The average (. . .) x can be calculated using the 
fact that the probability distribution P(x(t),t) of the Wiener process from Eq. (146]) is 
Gaussian: 

P( x ,t) = -=L=e-^£. (49) 
V 2ira z t 

It should be evaluated for a path x at two positions t and t' . Let us focus on t' > t, the 
case t' < t can be obtained by exchanging t ■<-> t' , n ■<-> m. The path to x(t') can be split 
into two parts, the path to x(t) and the path from x(t) to x(t'), where the difference 
x(t') — x(t) is a new random variable y(t' — t) that follows a Wiener process (146j) . but 
during the period if — t. This possible splitting is an expression of the fact that a Wiener 
process has no memory to the history of how the system reached the point x(t). We 
have 

iLJo[nx{t)+mx(t')}\ / iu)o[(n+m)x(t)+my(t' — t)] 



e 

\ I x,y 

e iu) {n+m)x(t)\^ ^iu my(t'-t) ^ 

where y{if — t) has the probability distribution 

1 (v-(t'-t)f 

*'-*)= /o fff, (51) 
^2-Ka z (t' — t) 
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The averages in Eq. (I50|) can be then written as Gaussian integrals: 



/ e iu (n+m)x(t)\ / e iw my(t' -t)\ _ 1 . / e ~ %Jl +ia) o (n+m)x ^ 
\ /*\ /w V27TO- 2 t 7-00 



oo 

2 



1 Z" 00 (y-(i'-t)) 

1 / e-W^- 4 *^*^, (52) 

oo 



^27rcr 2 (t'-t) 

and they can be easily calculated. Inserting the result back into Eqs. (148jl and (1471) 
and performing the integrals over t, t' we obtain (after some tedious calculations) the 
spectrum ( 1391) . 
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